Thermalization in one- plus two-body ensembles for dense interacting boson systems 
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Employing one plus two-body random matrix ensembles for bosons, temperature and entropy are 
calculated, using different definitions, as a function of the two-body interaction strength A for a 
system with 10 bosons (m = 10) in five single particle levels (N = 5). It is found that in a region 
A ~ At, different definitions give essentially same values for temperature and entropy, thus defining 
a thermalization region. Also, (m, N) dependence of A t has been derived. It is seen that At is much 
larger than the A values where level fluctuations change from Poisson to GOE and strength functions 
change from Breit-Wigner to Gaussian. 
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I. INTRODUCTION 



In recent years, the study of thermalization in iso- 
lated finite many-body quantum systems, due to inter- 
particle interactions, has received considerable interest 
[l|-l7|. This interest has arisen mainly due to major de- 
velopments in experimental study of many-particle quan- 
tum systems such as ultracold gases trapped in optical 
lattices 0. In particular, Rigol and Santos group em- 
ployed interacting spin- 1/2 systems (fermions and hard 
core bosons) on a lattice and examined various issues 
such as the role of localization and chaos, statistical re- 
laxation, eigenstate thermalization, ergodicity principle 
and so on [l|-|g|. Horoi et al. @ and Kota and Sahu [10| 
examined occupancies and different definitions of entropy 
for 28 Si and 24 Mg respectively, using nuclear shell model 
with realistic interactions. It was shown that the nuclei 
would be in the thermodynamic regime in general, when 
these are away from the ground state. Similarly, Casati's 
group examined in the past, Fcrmi-Dirac (FD) repre- 
sentability of occupation numbers for a spin- 1/2 system 
on a two-dimensional lattice 11] and more recently stud- 
ied thermalization in spin-1/2 systems, locally coupled to 
an external bath, using an approach based on the time- 
dependent density-matrix renormalization group method 
12]. Let us add that, as emphasized by Santos et al. @, 
it is possible to realize interacting spin-1/2 models ex- 
perimentally in optical lattices. 



On the other hand, thermalization in fcrmionic sys- 
tems has been studied in some detail using embedded 
Gaussian orthogonal ensemble of one- plus random two- 
body matrix ensembles [called EGOE(l+2)] [HEi], i.e. 
random matrix ensembles in many-fermion spaces gen- 
erated by random two-body interactions in presence of 
a mean-field. The EGOEs form generic models for fi- 
nite isolated interacting many-fermion systems (for bo- 
son systems these are called BEGOE with 'B' for bosons) 
and they model what one may call quantum many-body 
chaos [lj, [HI • The role of interactions in thermalization 
can be investigated by varying the interaction strength in 
these models. For example, for spin-less fermion systems, 



Flambaum et al. showed that EGOE(l+2) exhibits a re- 
gion of thermalization [l5[ and found the criterion for the 
occupancies to follow FD distribution jl5|, 16]. Later, us- 
ing EGOE both for spin-less fermions and fermions with 
spin [l(| Gjl Gil ; thermalization region generated by ran- 
dom interactions has been established by analyzing dif- 
ferent definitions for entropy. Going beyond these, in a 
more detail study, thermalization has been investigated 
within EGOE(l+2) for spin-less fermions by Kota et al. 
using the ergodicity principle for the expectation val- 
ues of different types of operators. Recently, Santos et al. 
compared results of spin-less EGOE(l+2) for statistical 
relaxation 19] with those from spin-1/2 lattice models 

m . 

Turning to interacting boson systems, thermalization 
was investigated by Borgonovi et al. [20] using a simple 
symmetrized coupled two-rotor model. They explored 
different definitions of temperature and compared the 
occupancy number distribution with the Bose-Einstein 
(BE) distribution. They conclude that: "For chaotic 
eigenstates, the distribution of occupation numbers can 
be approximately described by the BE distribution, al- 
though the system is isolated and consists of two particles 
only. In this case a strong enough interaction plays the 
role of a heat bath, thus leading to thermalization" . As 
BEGOEs [2l|-[25j are generic models for finite isolated 
interacting many-boson systems, it is important to in- 
vestigate thermalization using these ensembles. 

Embedded Gaussian orthogonal ensemble of one- plus 
random two-body matrix ensembles for spin-less boson 
systems is called BEGOE(l+2) and this ensemble was 
introduced and analyzed for spectral and wave-function 
properties in |2l| - [23| . For m bosons in N single par- 
ticle (sp) levels, in addition to dilute limit (defined by 
m — y 00, N — > 00 and m/N —} 0), another limiting situa- 
tion, namely the dense limit (defined by m — > 00, N — > 00 
and m/N — > 00) is also feasible. This limiting situation 
is absent for fermion systems. Therefore the focus was on 
the dense limit in BEGOE investigations [2l| - [25j . In the 
strong interaction limit, two-body part of the interaction 
dominates over one-body part and hence BEGOE(l+2) 
reduces to BEGOE(2). Some of the generic results es- 
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tablished for BEGOE(l+2) are as follows: (i) eigenvalue 
density approaches Gaussian form [2l|, [26[ ; (ii) for strong 
enough interaction, there is average-fluctuation separa- 
tion in eigenvalues [2l|, [I?) ; (hi) similarly, the ensemble 
is ergodic in the dense limit with sufficiently large N 
[22] | and there will be deviations for small N (25|; (iv) 
as the strength of the two-body interaction, A, increases, 
there is Poisson to GOE transition in level fluctuations at 
A = A c [22j and with further increase in A, there is Breit- 
Wigner to Gaussian transition in strength functions at 
A = Xp [HI ■ The main result of the present Letter is the 
demonstration that finite dense interacting boson sys- 
tems generate a third chaos marker At [as in fermionic 
EGOE(l+2) ensembles], a point or a region where differ- 
ent definitions of entropy, temperature, specific heat and 
other thermodynamic variables give the same results, i.e. 
where thermalization occurs. 



II. TEMPERATURE: DEFINITIONS AND 
RESULTS 

Temperature can be defined in a number of different 
ways in the standard thermodynamical treatment. These 
definitions of temperature are known to give same result 
in the thermodynamical limit i.e. near a region where 
thermalization occurs [l|. In this section, four differ- 
ent definitions of temperature (T = described be- 
low, have been used to compute the temperature of finite 
dense interacting boson systems as a function of energy 
as well as a function of the two-body interaction strength. 

• j3 c : defined using the canonical expression, be- 
tween energy and temperature which allows stan- 
dard thermodynamical description for the quantum 
system, is given by 



In this Letter, we present results for thermalization in 
dense interacting bosonic systems by varying the strength 
parameter A of the two-body interaction in the BE- 
GOE(I+2) Hamiltonian given by H = h(l) + XV (2). 
Here h(l) is one-body part of the interaction, defined 
by single particle energies (SPEs) £k (k = 1 to N) for N 
sp levels while the two-body interaction V(2) is defined 
by the two-body matrix elements (TBMEs) denoted as 
Vi jM = {(ij)\V{2)\(kl)}. In the present study SPEs are 
taken as independent gaussian random variables with 
mean equal to k and variance equal to 1/2. Similarly, 
TBMEs are taken as independent gaussian random vari- 
ables with zero mean and variance=l for off diagonal 
TBMEs and variance=2 for diagonal TBMEs. Construc- 
tion of the m-boson Hamiltonian, H(m), and thereby 
the BEGOE(l+2) ensemble in m-particle space with ma- 
trix dimension d = ( N+ ™~ 1 ) was described completely in 
[2l|,[23j]. The results, presented here, have been obtained 
by fully diagonalizing 100-members of a BEGOE(l+2) 
ensemble with 10 bosons in 5 sp levels for each value of A. 
The dimensionality of the system is d = 1001. (we have 
also carried out calculations for 10 bosons in 4 sp levels 
and similar results were obtained, but they are not pre- 
sented here as this is a much smaller example, d = 286). 
The ensemble average is carried out by making the spec- 
tra of each member of the ensemble zero centered (e is 
centroid) and scaled to unit width (a is width). 
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where Ei are the eigen-energies of the Hamiltonian. 
With above relation, (E) can be obtained at given 
j3 c using all the eigen-energies of the system. 

(5 fa: defined using occupation numbers obtained 
by making use of the standard canonical distribu- 
tion is given by, 



J2i n ( E i)k exp[-/3/j t Ei 
Ei exp[-0/« Ei] 



(2) 



Here k is sp level index and i is eigen-energy index. 
Using expectation values of occupancies calculated 
for all eigen-states and exact eigen-energies, j3fa 
can be computed by considering Eq. (2) as one pa- 
rameter fitting expression and with the constraint, 

• Pbe'- defined using BE distribution for the occupa- 
tion numbers is given by, 



n(£)£ E = l/{exp[pBE(E) (e k - /*(£))] - 1}. 



(3) 



The paper is organized as follows. In Section HH differ- 
ent definitions of temperature are given and results ob- 
tained by varying the two-body interaction strength in 
BEGOE(l+2) are described. Similarly, results obtained 
using three different definitions for entropy, are described 
in Section Hill They allow us to define the thermalization 
marker At. In Section HVl duality point is discussed, us- 
ing information entropy and strength functions, in the 
two extreme basis defined by h(l) and V(2) operators 
and the (m, N) dependence of the marker At is derived. 
Finally, Section [V] gives conclusions. 



Here fj, is a chemical potential. Although, this ex- 
pression is derived for many-body non-interacting 
particles in contact with a thermostat, it is shown 
that conventional quantum statistics can appear 
even in isolated systems with relatively few par- 
ticles, provided a proper renormalization of en- 
ergy is taken [l5|, [lg] . Comparing numerical data 
of expectation values of occupancies at a partic- 
ular eigenenergy and given SPEs, the unknowns 
Pbe and H in the BE distribution with constraint 
Efc n (E) BE — m can be obtained. 
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• Pt- defined using state density, p(E), of the total 
Hamiltonian. The thermodynamic entropy is de- 
fined as S ther {E) = \og\p{E)]. The j3 T can be com- 
puted using Pt = dXn ]iE E ^ ■ Here not that, for BE- 
GOE(l+2), the form of state density is very close 
to Gaussian irrespective of the value of two-body 
interaction strength [2l], [26[ . 

Figure [T] shows ensemble averaged values of /?, com- 
puted via definitions described above, for a 100 member 
BEGOE(l+2) ensemble with m = 10 and N = 5 as a 
function of normalized energy, E — (E — e)/a, for var- 
ious A values. Here, we compare numerical values of p 
from E = — 1.5 to the center of the spectrum, where tem- 
perature is infinity. The edges of the spectrum have been 
avoided for the following reasons: (i) density of states is 
small near the edges of the spectrum; (ii) eigenstates near 
edges are not fully chaotic. The j3 values are obtained 
for all members separately and then ensemble average is 
carried out taking bin-size equal to 0.1. Since the state 
density for BEGOE(l+2) is Gaussian irrespective of A 
values, Pt a s a function of energy gives straight line. 
In Fig. [TJ Pt results are shown in the plots by dotted 
lines. It is clearly seen from Fig. Q] that for the inter- 
action strength A < A c (For (m, N) = (10,5), A c ~ 0.02 
and Xp ~ 0.05, see ref. 22]), there is significant differ- 
ence between the numerical values of P obtained via var- 
ious definitions of temperature. Going further beyond A c 
(A c < A < Xp) where GOE fluctuations in state density 
sets in but the eigenstates are still not fully chaotic, the 
P values obtained via canonical expressions, defined by 
Eq.(l) and (2), give good agreement near the center of 
the spectrum. There are deviations near low temperature 
region. Near the region X — Xp and beyond, the eigen- 
states become fully chaotic giving very good agreement 
between the numerical values of p c and Pfu- The inverse 
temperature Pbe is obtained by solving microcanonical 
definition given by Eq. (3), with SPEs taken as indepen- 
dent Gaussian random variables and results are shown 
in the plots by red stars. In the region X < Xf, inverse 
temperature Pbe, found from BE distribution turns out 
to be completely different from /3 values obtained using 
other definitions. As in this region, the structure of eigen- 
states is not chaotic enough, leading to strong variation 
in the distribution of the occupation numbers and thus 
strong fluctuations in Pbe- Moreover, near the center 
of the spectrum (i.e. as T — > oo), the value of denom- 
inator in Eq. (3) becomes very small, which leads to 
large variation in Pbe values from member to member. 
Further increase in A > Xp, in the classically chaotic re- 
gion, the occupation number distribution becomes statis- 
tically stable with respect to the choice of eigenstate and 
at one point A = At, temperatures defined using canon- 
ical and microcanonical definitions give same result, i.e. 
Pfit ~ Pbe- This lead to same values of temperature 
giving the thermodynamic marker At. In the strong in- 
teraction domain (A > At), the match between Pbe with 
other P values is not good. This is due to neglect of 



induced SPEs, = J2j Vkjkj , from two-body inter- 
action [26]. When the interaction is weak, induced SPEs 
part is small but in strong interaction domain their con- 
tribution is important. Adding induced SPEs part into 
SPEs, Pbe is obtained by taking a proper renormaliza- 
tion of energy and results are also presented in Fig. 1 
for A values 0.13 and 0.2 by green filled circles. Here, 
Pbe values come quite close to other p values. We found 
that the match between different values of P is good near 
X = Xt — 0.13. In the next section, we present our re- 
sults for similar study using the different definitions of 
entropy. 



BEGOE(l+2) : m=10, N= 5 
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FIG. 1. (Color online) Ensemble averaged values of inverse 
temperature (/3) as a function of normalized energy, (E — 
e)/cr, for different values of two body interaction strength A, 
calculated using a 100 member BEGOE(l+2) ensemble with 
m = 10 and TV = 5. Here in calculations, different definitions 
of inverse of temperature are used. In the calculations, SPEs, 
drawn from independent real Gaussian random variables, are 
employed. Results are also shown, in the plots for A = 0.13 
and 0.2, for Pbe by green filled circles, in which induced SPEs 
from the two-body interaction are taken into account. 
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III. ENTROPY: DEFINITIONS AND RESULTS 

In this section, to identify the thermalization region 
and hence value of the third marker At for BEGOE(l+2), 
we consider following three definitions of entropy. 

• Thermodynamic entropy, obtained using the state 
density p{E), as a function of of energy eigenvalues; 
S ther {E) = log[p(£)]- 

• Information entropy in the mean-held basis defined 
by S m f°(E) = - Y?i\Cf | 2 log(|Cf | 2 ), here |Cf | 2 
is the probability of basis state i in the eigenstate 
at energy E. 

• Single particle entropy, obtained by calculating 
the occupancy of different single-particle states, 
as a function of energy eigenvalues; S sp (E) — 
~ J2k ( n k{E)) \og((rik(E))). Here the summation 
is over all N sp levels and nu{E) is the occupancy 
of the fc-th sp level at energy E. 

We use following measure, defined using above defini- 
tions of entropy [7j, to obtain A*: 



A.(A) = 



R thery 



(R s i 



Rf er ) 2 ]dE 



JZoRT r dE 

(4) 

where R% = exp[S a (E) - S* ax ]. In the thermody- 
namic region the values of the different entropies will be 
very close to each other, hence the minimum of A s gives 
the value of X t . In Fig. [21 results shown for ensemble 
averages A S (A) (blue stars) obtained for a 100 member 
BEGOE(l+2) ensemble with 10 bosons in 5 sp levels as 
a function of A. The second vertical dash-line indicates 
the position of At where ensemble average A S (A) is min- 
imum. For the present example, we obtained A t ~ 0.13. 
This value of At is same as obtained in Section|TTl at which 
different definitions of temperature give same values. 

In the past it is demonstrated that for BEGOE(l+2), 
as the strength A of the two-body interaction increases, 
there is Poisson to GOE transition in level fluctua- 
tions at A = A c [llj]. In order to show that A c << 
A t , we study the nearest neighbor spacing distribution 
(NNSD) as a function of A to detect the position of 
marker A c . It is well known that when the system 
is in integrable domain, the form of NNSD is close to 
the Poisson distribution, i.e, P(s) = exp(— s), while in 
chaotic domain, the form of NNSD is the Wigner sur- 
mise, P(s) = {tts/2) exp(— 7rs 2 /4). To interpolate be- 
tween these two extremes, we use Brody distribution f28j j . 
P{s,oj) = A u (uj + l)s w exp(-A; s" +1 ). Here w is called 
the Brody parameter and A u is a normalization constant. 
If the spectral fluctuations are close to the Poisson type, 
u = 0, or to the Wigner surmise with cj = 1. The po- 
sition of the chaos marker A c is fixed by the condition 
w(A) = 1/2. Here, the NNSD is obtained, using the 



unfolding procedure described in [2l|, with the smooth 
density taken as a corrected Gaussian with corrections 
involving up to 6th order moments of the den sity fu nc- 
tion. In Fig. [21 Ensemble averaged values of w(A) are 
shown by filled triangles for a 100 member BEGOE(l+2) 
ensemble with (m, jV)=(10,5) as a function of two-body 
interaction strength A. The A c is shown by vertical dash 
line in the Fig. [21 Here, different criterion is used, than 
in [22J, to obtain A c although match is very good. The 
results clearly show that A c << A t for BEGOE(l+2) just 
as seen for EGOE(l+2) fermionic ensembles 0, GI- 
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FIG. 2. (Color online) Ensemble averaged values uJ and A 3 
as a function of two body interaction strength A, calculated 
using a 100 member BEGOE(l+2) ensemble with m — 10 
and N = 5. The vertical dash-lines represent the position of 
A c and At. Here A c ~ 0.02 and At ~ 0.13. 



IV. DUALITY AND (m, N) DEPENDENCE OF A t 

The duality region A = A^ is the point where quanti- 
ties defining eigenstate properties like entropy, strength 
functions, temperature etc. give same values irrespec- 
tive of the defining basis [29|]. Then one can argue 
that in this region all eigenstates look alike and the 
duality region defined by A — Xd is expected to corre- 
spond to the thermodynamic region defined by A = At 
as in fermionic EGOE(l+2) results [13, EI- For the 
BEGOE(l+2) Hamiltonian, two choices of basis appear 
naturally. One is mean-field basis defined by h(l) and 
another is the infinite interaction strength basis defined 
by V(2). To examine duality, we compare information 
entropy S m f°(E) and strength functions F^ k (E) (also 
called local density of states (LDOS)) in h(l) basis and 
in V(2) basis. The strength function corresponding to 
the fc'th basis state for a particular member of ensem- 
ble is defined by F ik (E) = J2i\ C k\ 25 ( E ~ E i), where 
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fc-energies = (k\H\k) and is the £;'th basis state 
for m-particles in N sp states. Figure [3] shows numerical 
results for S info {E) and F ik (E) in h(l) and in V{2) basis 
for a 100 member BEGOE(l+2) ensemble with (to = 10, 
N = 5) for different values of A. Here strength functions 
Ff u (E ) are computed following the procedure described 
in [23[ and results are shown for £/. = in the Fig. \5jp. It 
is seen from Fig. [3] that values of S mfo (E) and % (E 1 ) in 
these two basis are found very close near A = 0.13 giving 
value for the duality marker A^ ~ 0.13 for the present ex- 
ample. For A < X d , the S lnfo {E) values in the h(l) basis 
are smaller compared to those in the V(2) basis and for 
A > Xd, S ln f° in the ft-(l) basis is comparatively larger. 
While opposite behavior is observed from the results of 
the strength functions. The values of entropy as well 
as of strength functions in these two basis coincide near 
A = Xd- This value is very close to the marker At and 
therefore, Xd region can be interpreted as the thermody- 
namic region in the sense that all different definitions of 
temperature and entropy coincide in this region. 

In the h(l) basis, S info (E) is determined using equa- 
tion given below [22| • 



expfS" 



info 



(E)-S%&] = exp{^} exp{-^} 

(5) 

Here £ is the correlation coefficient between the full 
Hamiltonian H and the diagonal part of the full Hamil- 
tonian H ; it is given by 



a 



o f f — diagonal 



(J 



H 



\ 



'h(l) 



+ A2 4( 2) 



(6) 



'Mi 



We can determine the value of At by using the con- 
dition that £ 2 = 0.5 [171; i.e. the spreadings produced 
by h(l) and V(2) are equal at At- In Fig. [4] ensem- 
ble averaged values of C 2 as a function of A for a 100 
member BEGOE(l+2) ensemble using (to, 2V)=(10, 5) is 
presented by filled red circles. It is clear from the fig- 
ure that for A < A c , £ 2 is close to 1 and as A increases, 
£ 2 goes on decreasing smoothly. The two vertical dash- 
lines in Fig. |31 indicate the respective positions of A c 
and At as obtained in Section IIII1 It can be clearly seen 
that ( 2 = 0.5 gives the thermalization point A f = 0.13. 
For BEGOE(l+2) ensemble, analytical expression for ( 
based on the method of trace propagation is derived in 
[13. With C 2 = 0.5 in Eq.(7) of [23] and solving it for A, 
(to, N) dependence of marker is given by, 



(N + 2)X 



N(N + l)(N - 2)(to - 1)(N + to + 1) 



(7) 



where X 
with Si = i, 



X 



E 



N -2 

to . 



N(N 



For uniform sp spectrum 
+- 1)(N - 1)/12 and At = 



(jV-l)(jV+2) 



3(JV-2)(m-l)(JV+m+ 



jj. With (to = 10, N = 5), we have 



At ~ 0.15. For single particle energies which we have 
used in the present study, X = N(N 2 + 5)/12 and 



X( 



(N + 2)(N 2 + 5) 



3(N + 1)(N - 2)(m - 1)(N + to + 1) 



(8) 



With (to = 10, N — 5), we have A t « 0.16. Figure 
4 shows plots of ( 2 as a function of A obtained using 
Eq.(7) of [13] ■ The blue curve in the figure is obtained 
due to uniform SPEs while the green curve is obtained 
due to SPEs employed in the present study. It can be seen 
from results that the ensemble averaged values are close 
to the expected values. Small discrepancy is due to the 
neglect of induced single-particle energies. In the dense 

limit, Eq.© gives A t ~ m\p^- Similarly, in the dilute 

limit, we have At ~ ,„ and this result is in agreement 

V3m , , ° 

with EGOE(l+2) result given in [l7j. From Eq.© it is 
seen that for m/N fixed as to — » oo and N oo (also 
into strict dense limit), A* — >■ 0. A similar behavior is 
expected for A c and Xf- These sudden transitions with 
A > are similar to the situation with Poisson to GOE 
or GUE HEll and GOE to GUE [H. 



O K(2) 




FIG. 3. (Color online) Ensembles averaged results for (a) 
the information entropy and (b) strength functions F^ k (E) 
in the h(l) and V(2) basis for a 100 member BEGOE(l+2) 
ensemble with (m = 10) and (N = 5) are shown as a function 
of normalized energy, (E—e)/a, for different A values. Results 
averaged over bin-size 0.1 are shown as circles; filled circles 
correspond to h(l) basis and open circles correspond to V(2) 
basis. Ensemble averaged £ values are also given in the figure. 
Strength function plots are obtained for = and in the 
plots / F^ k (E)dE = 1. 



6 



1.0 



X 

c 


V\ X 

Va t 

^ 














■ -•-BEGOE(l+2) 




Fix SPE 




Random SPE 





1E-3 0.01 0.1 1 



X 

FIG. 4. (Color online) Ensemble averaged values of £ 2 (red 
filled circles) as function of the two-body interaction strength 
A, calculated for a 100 member BEGOE(l+2) ensemble with 
(m = 10, N — 5) are shown. Smooth curves are obtained 
using Eq.(7) of ref . [23| . Blue (continuous) curve is for fixed 
SPEs and green (dash) curve is due to random SPEs used in 
present study. 



V. CONCLUSIONS 

In the present work, we have analyzed the relation- 
ship between order to chaos transition and thermaliza- 



tion in finite dense interacting boson systems using one- 
plus two-body embedded Gaussian orthogonal ensemble 
of random matrices. Using numerical calculations, it is 
demonstrated that in a region A ~ A t different defini- 
tions give essentially same values for temperature and 
entropy, thus defining a thermalization region similar to 
as in fermionic EGOE(l+2) ensembles. The value of At 
is much larger than the A value where level fluctuations 
change from Poisson to GOE and strength functions from 
Breit-Wigner to Gaussian (A t > Xf > A c ). Further it is 
established that the duality region where information en- 
tropy and strength functions will be the same in both the 
mean field and interaction defined basis corresponds to a 
region of thermalization. We have also obtained formula 
for At in terms of (m,N). In addition to this, we know 
that At >> A c , which is further confirmed by the analyt- 
ical formula for At given in Eq.© and the estimates for 
A c as given in [22j . Similar results are known for fermion 
where A c oc l AT and At oc -j=. However, for bosons 
formulas are not available for A c and A^. Therefore we 
cannot tell if A^ will be close to A c or A t or it will be 
far from both. This is an important open problem. The 
present work brings completion to the study of transition 
(chaos) markers generated by BEGOE(l+2) initiated in 
[23, l23| . Further investigations on thermalization in BE- 
GOE(l+2) will be discussed in future. 
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